Charge and spin transport in strongly correlated one-dimensional quantum systems 

driven far from equilibrium 
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We study the charge conductivity in one-dimensional prototype models of interacting particles, 
such as the Hubbard and the t — V spinless fermion model, when coupled to some external baths 
injecting and extracting particles at the boundaries. We show that, if these systems are driven 
far from equilibrium, a negative differential conductivity regime can arise. The above electronic 
models can be mapped into Heisenberg-like spin ladders coupled to two magnetic baths, so that 
charge transport mechanisms are explained in terms of quantum spin transport. The negative 
differential conductivity is due to oppositely polarized ferromagnetic domains which arise at the 
edges of the chain, and therefore inhibit spin transport: we propose a qualitative understanding of 
the phenomenon by analyzing the localization of one-magnon excitations created at the borders of a 
ferromagnetic region. We also show that negative differential conductivity is stable against breaking 
of integrability. Numerical simulations of non-equilibrium time evolution have been performed by 
employing a Monte-Carlo wave function approach and a matrix product operator formalism. 

PACS numbers: 75.10.Pq, 05.30.-d, 05.60.-k, 03.65.Yz 



O 
o 



> 

CO 

O 

CN 

o 

On 
O 



X 



I. INTRODUCTION 

Transport properties of strongly interacting fermions 
in microscopic models of one dimensional quantum sys- 
tems have been the subject of a large number of theo- 
retical and experimental studies^. In the last few years 
this has become a topical subject, due to the rapidly 
developing process of miniaturization in semiconductor 
microelectronics devices that is approaching its natu- 
ral limits, reaching the atomic or molecular scale^^. 
Of course, a technological breakthrough in this direc- 
tion would require conceptually new devices, such as few 
or even single molecules embedded between electrodes, 
that could perform the basic functions of microelectron- 
ics. The first promising step in the realization of such 
devices comes from the observation of many-body ef- 
fects, as the Coulomb blockade and the Kondo effect in 
nanometer-scale systems, like single molecules or carbon 
nanotubes 5 -^. Establishing and reaching a suitable de- 
gree of control of nonlinear electronic transport, such 
as a Negative Differential Conductivity (NDC) regime, 
would be one of the ultimate tasks for functional nan- 
odevices, since it lays at the basis of current rectification 
and amplification. NDC has been observed in a variety 
of nanoscopic objects, like semiconductor quantum dots 7 , 
carbon nanotubes 8 -, as well as single molecules^. 

From a theoretical point of view, nonlinear transport 
properties in such systems are usually studied by con- 
sidering effective models of few single-particle levels, see 
e.g. Refs. [lOliTllfTifTlili. In this paper we adopt a rather 



different perspective and study the full many-body quan- 
tum dynamics of one-dimensional prototype models of 
strongly interacting fermions, when they are coupled to 
some external baths. We will show how effects of non- 
linear transport naturally emerge in far-from-equilibrium 
situations, by exploiting the many-body dynamics of such 
microscopic models in its whole complexity. While sit- 
uations close to equilibrium are quite well understood 
and can be tackled by the powerful linear response for- 
malis m 15 i 16 i 17 i 18 ' 19 i 20 ' 21 , almost nothing is known about 
the physics of such systems far from equilibrium. In 
this regime new quantum phases and phenomena can ap- 
pear, thus making the problem relevant also for funda- 
mental physics 22 . Furthermore, the study of far-from- 
equilibrium quantum systems is of interest also for issues 
such as the control of heat flow at the nanoscal o 23 ' 24 and, 
in quantum information processing, for quantum state 
preparation/transfer 25 . Unfortunately a fully analyti- 
cal treatment is generally unfeasible 26 , and one typically 
has to resort to numerical simulations, aimed at solving 
the quantum master equatio n 27 ' 28 i 29 ' 30 i 31 ' 32 i 33 , or based 
on different approaches, like path integral Monte Carlo 
approach 34 , time-dependent density matrix renormaliza- 
tion group or current density functional theor y 35 i 36 i 37 , 

In this paper we consider two prototype microscopic 
one-dimensional models of interacting fermions, namely 
the Hubbard model and the t — V spinless fermion model, 
and couple them to some external baths that inject and 
extract particles at the system edges, thus mimicking the 
effect of electrodes. The Hubbard and the t — V model 
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can be mapped into the Heisenberg spin- 1/2 ladder and 
chain, respectively. In these spin models NDC reflects in 
the suppression of spin conduction, while the operators 
injecting/extracting electrons are mapped into operators 
flipping the two spin species at the border of the chain. 
As an example of spin chains coupled to such "magnetic 
baths" one can consider molecular spin wire s 38 ! 39 with 
each boundary coupled to an external spin (magnetic 
impurity); the ratio of up/down and down/up spin- flip 
probabilities is determined by the populations of such im- 
purities, which in turn can be tuned by means of applied 
electromagnetic fields. In the linear response regime, the 
electronic (i.e., fermionic) transport and correspondingly 
the spin transport can be ballistic or diffusive, depend- 
ing on the values of the Hamiltonian parameters. Here 
we focus on the far- from- equilibrium regime, beyond the 
linear response regime. Our numerical results show that, 
strikingly, in the above mentioned models it is possible 
to achieve a regime where charge/spin conductivity ex- 
hibits a negative differential with respect to the driving 
strength. NDC arises as a result of the appearance of 
a far-from-equilibrium steady state characterized, for the 
spin chain models, by long-range spin ordering into ferro- 
magnetic domains. These ferromagnetic domains corre- 
spond to charge separation in the fermionic models, with 
all the electrons frozen in half of the lattice. In both 
cases, it is clear that such cooperative many-body state 
hampers spin flips (or charge injection/extraction), thus 
strongly suppressing the current. We will show that our 
numerical results can be qualitatively explained in terms 
of localization of one-magnon excitations. 

The paper is organized as follows: in Sec. [TT] we start 
by setting our electron transport problem and reducing 
it to a Lindblad master equation formalism, that will be 
used throughout the paper. In Sec. IIIII we introduce the 
model of open Hubbard chain coupled to two macroscopic 
reservoirs, and discuss some peculiar charge transport 
properties, focusing on the NDC behavior. In Sec. IIVI 
we consider a simplified model for spinless fermions and 
show that it can be mapped into a Heisenberg spin chain. 
In Sec. [V] we study in details the spin transport prop- 
erties of the Heisenberg chain. Moreover, we provide 
a one-magnon localization argument which qualitatively 
explains the observed NDC behavior. To explore the 
possibility that our system undergoes a metal-insulator 
phase transition when driven far from equilibrium, we 
propose to study steady-state spin-spin correlation func- 
tions. We also add, in Sec. [VI] a staggered magnetic field 
and check that NDC is stable against breaking of integra- 
bility. Finally, in Sec. ELD we draw our conclusions. In 
the Appendices we describe the two numerical methods 
used throughout the paper, namely the quantum trajec- 
tories approach and the matrix product operator formal- 
ism (App. [A"]l . give technical details on the mapping of 
our fermionic systems into spin chain models (App. |B|) . 
provide some numerical results about the steady-state 
spin-spin correlation functions fApp. [C]). and present an 
analytical derivation of the one-magnon argument for the 




FIG. 1: Schematic drawing of the level structure for a chain 
with N — 3 sites. We assume that the lead-chain tunneling 
rates Fl,!^ are much larger than the intra-chain tunneling 
rates fii, f^- The electronic current flows from the right lead 
(emitter) to the left lead (collector). 



Heisenberg spin chain (App. [DJ. A brief account of the 
NDC features of the Heisenberg chain can be found in a 
recent paper by some of us 3 ^. 

II. MASTER EQUATION APPROACH 

Our electronic transport model is described by the 
Hamiltonian 

H = Hs + Hi+H c , (1) 

where the different terms correspond to the nanoscale 
electronic system, the leads, and the leads-system cou- 
pling, respectively. 

As sketched in Fig. [I] we consider a iV-site chain, whose 
autonomous dynamics is described by Hamiltonian tis- 
Such a lattice models a nanoscale system, for instance a 
chain of coupled quantum dots or a molecular wire (in the 
latter case, each lattice site corresponds to one atom). 

The first and the last site of the chain, 1 and N, are 
coupled to the left /right leads via the tunneling Hamil- 
tonian 

H c = J2(T L kc{ Ks c hs + T Rk c Bks c N , s ) + H.c, (2) 

fc,s 

where c\ c are fermionic creation/annihilation operators: 
cj s creates an electron with spin s at site j (j — 1, N, 

s =t)4)) c Lfc s ( c flfcs) crea t es an electron in the left 
(right) lead in the state \Lk, s) {\Rk, s)). 
The leads are modeled as ideal Fermi gases, 

H i = e fc( c L,s Cifc > s + 4 M c flM ), (3) 

k,s 

which arc initially at equilibrium, at temperature T and 
chemical potentials /il, [Mr = fir, + eV, where V is the 
applied bias voltage and e the electron charge. We as- 
sume that the coupling between the system and the leads 
is weak, such that the state ps(t) of the leads at any time 
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t is well described by pB (t) = pl ® PR, with p^ and p R 
grand canonical density matrices for the left and right 
leads, respectively. 

A Lindblad master equation for the system's evolu- 
tion can be obtained from our microscopic Hamiltonian 
model following standard textbook derivations, under the 
usual Born-Markov and rotating wave approximations 
and neglecting the Lamb-type renormalization of the un- 
perturbed energy levels (see, e.g., Sec. 3.3 of Ref. liol): 



dp 

at 



:[U S ,P] 



.P L tn 



(4) 



where p(t) is the density matrix describing the open 
quantum system, the Lindblad operators L m describe 
the effect of the environment, while [•, •] and {•, •} denote 
the commutator and the anti-commutator, respectively. 
Hereafter we shall set H = 1. Moreover, we assume that 
the tunneling between sites 1 (N) and the left (right) 
lead is much faster than the intra-chain tunneling and 
that we can neglect the effects of Coulomb repulsion on 
the system-leads transition rates. This is the case in the 
so-called wide-band limit, in which the conduction band- 
width of the leads is much larger than all other relevant 
energy scales and all the relevant lead states are located 
in the center of the conduction band, so that the energy 
dependence in the system-leads transition rate may be 
neglected. Under these approximations, we can easily 
derive the Lindblad master equation ([4]) , with four Lind- 
blad operators on each of the two chain ends: 



L 3 = VTlJZc f u 
and similarly 
i 5 = vTXfec t 



Li 
L 4 



h) ci.t , 



jv,t 

L7 = \Z^rJr, cjv_| 
where 



T L = 2n\T L (E 1 )\ 3 g L {E 1 ), f L = f L {E x ) 



(5) 



(6) 



(7) 



T R = 2ir\T R (E N )\ 2 g R (E N ), f R = f R (E N ), (8) 

with Ei (En) being the energy difference between the 
two chain states involved in the transitions for site 1 
(N), gi (I = L,R) being the density of states of lead I 
(we assume that the leads are macroscopic objects, with 
a continuous density of states), Tj(e = efc) = T/fe, and 

//(e) = [l + e( e ~ Mi '/ fcBT ] denoting the Fermi function, 
with ks the Boltzmann constant. Note that the energy 
differences E\, En contain the charging energy E c , if an 
electron is tunneling onto an already occupied site, but 
does not contain it if the site is initially empty. We have 
neglected the dependence of the Fermi functions 
on E c . This condition is fulfilled when E c <C k R T (in the 
Hubbard and t — V models described in this paper such 



constraint corresponds to on-site repulsion U <C fcsT 
and nearest-neighbor repulsion V <C fcgT, respectively). 
Finally, in order to consider incoherent tunneling of elec- 
trons into the chain (sequential tunneling approxima- 
tion), the level broadening due to the chain-leads tun- 
neling must be small compared to temperature, that is, 
we require r^T^ < i; B T. 

As we shall see in this paper, a main advantage of the 
master equation approach is that it can be applied far- 
from-equilibrium, beyond linear response regime. The 
far-from-equilibrium regime in our model corresponds to 
large bias voltage, eV ksT, with the energy differ- 
ences Ei, En such that p^ <C E\,En <C p R - In this 
limit, f]j — > 0, f R — * 1, that is, the backward flow of 
electrons (against the applied bias) vanishes. 

The master equation approach may be generalized, in- 
cluding the effects of Coulomb repulsion^ (thus describ- 
ing the Coulomb blockade phenomenon) or the coupling 
of multilevel nanoscale systems to external lead o 42 ' 43 . 
The price to pay for such generalizations is, in general, 
the introduction of a larger number of Lindblad oper- 
ators, corresponding to all possible transitions between 
the levels of the nanosystem. 

The main transport quantity, the electron current j, 
is defined by the continuity equation of the local charge 



density rife , 



.Cfe )S , n k = n k ^ + nk,i- 

^ r + v, fe = o, 



which can be rewritten as 



- jk = i[nk,Us], k = l, 



,N-1. 



(9) 



(10) 



Note that, due to the continuity equation, one has j = jk 
for any k along the chain. 

By definition the electron current is given by 



(i) = 



dN R dN L 



dt 



dt 



(11) 



with N R (Nl) being the number of electrons in the right 
(left) lead. This equation expresses the current in terms 
of the number of electrons which enter the system from 
the left reservoir (—dN^/di), or go out into the right 
reservoir (dN R /dt) per time unit. As we shall discuss in 
Appendix El —dN^/dt and dN R /dt may be computed 
by means of the quantum trajectories approach. We will 
use both Eqs. (fTO]) and Eq. ifTT]) to compute the current. 

Explicitly solvable models of master equations are very 
limited, therefore support from extensive numerical sim- 
ulations is generally needed. We used two methods to 
face this problem. The first is a Monte-Carlo wave 
function approach, that is based on the technique of 
Quantum Trajectories (QT), widely used in quantum op- 
tic o 44 i 45 ' 46 i 47 . The second is a Matrix Product Operator 
(MPO) technique based on the time-dependent density 
matrix renormalization group metho d 48 i 49 ' 50 i 51 ' 52 . QT 
revealed themselves a powerful tool in the study of rela- 
tively small system sizes, especially for situations with a 
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strong external bias, where equilibration times needed to 
reach the steady state are generally long. On the other 
hand, the MPO method can deal with systems up to one 
order of magnitude larger, but it may encounter some 
difficulties in converging to the stationary state for large 
driving fields. Both numerical methods are briefly dis- 
cussed in Appendix [X] 



III. HUBBARD MODEL 

We start our analysis by considering a paradigmatic 
model for the physics of strongly interacting electronic 
systems: the Hubbard model. Its Hamiltonian is a sum of 
a kinetic term allowing for electron tunneling between the 
neighboring lattice sites, and a potential term consisting 
of an on-site interaction; in one dimension it is given by: 

N 

H S = -tJ2 (4,s c J+M +H.c.) +uY j n jA n jA , (12) 

2> 2=1 

where s stands for spin up f or down j, configuration, 
while j = 1, . . . , N is the site index and N is the number 
of lattice sites. The operators cj s , c JlS create/annihilate 
a spin- 1/2 fermion with spin s at site j, and satisfy the 
usual anticommutation rules; rij iS = ct s Cj >s is the cor- 
responding number operator. We consider open bound- 
ary conditions, therefore the sum over j in the first term 
runs from 1 to TV — 1. The system parameters t and U 
(U > 0) describe, respectively, the nearest neighbor hop- 
ping strength and the on-site repulsion between electrons 
with opposite spins. 

Both ends of the Hubbard chain are coupled to some 
electrodes which act on the system by injecting or ex- 
tracting particles with different spins. In the Lindblad 
master equation formalism, we assume that their effect 
can be modeled by the Lindblad operators ([5]) and ©■ 

From the continuity equation (| 10|) we obtain 

.? = -*I> c Lc fc+M +H.c.), (fc = l,...,JV-l). (13) 



We examined the fermionic transport properties of the 
Hubbard model (|12p coupled to external baths by exploit- 
ing a mapping of this system into a spin ladder model, 
where the particle current is replaced by the spin cur- 
rent. Specifically, the Hamiltonian in Eq. (|12[) is mapped 
into a Heisenberg spin ladder by first employing a double 
Jordan- Wigner transformation of spin-up and spin-down 
fermions (separately) into two different species of hard 
core bosons. Then they are transformed into two species 
of spin-1/2 particles, which are described by the Pauli 
matrices cr? and rj* (a = x,y,z). Details are given in 
Appendix [Bj One finally arrives at the following spin 
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FIG. 2: (Color online). Spin current for the a (full curves 
and symbols) and the r species (dashed curves and empty 
symbols) of spin as a function of the driving strength for the 
Hamiltonian in Eq. (|14|l with [7 = 5. The system-bath cou- 
pling is set equal to F = 0.5; the simulation time (QT ap- 
proach) is T = 2 x 10 5 . Note that curves and symbols for the 
a and r species are nearly superimposed. For the Hubbard 
model we set t = 1 as the system's energy scale. 



ladder Hamiltonian for the autonomous system: 
JV-l 

H * = - 2 £ I^ "- 1 + + + T ^ l} 

2=1 

TT N 

2=1 

The Lindblad operators in Eqs. (O, © correspond, in 
the spin-1/2 picture, to operators flipping the two spin 
species at the borders of the chain. Apart from some 
phase factor which is uninflucnt for our purposes (see 
Appendix |SJ we have that ct ^ — > and Cj.j — > aj for 
spin- up particles, while cj . — > and Cjj — > tJ for spin- 
down particles [<jj = (cr* ±icxJ)/2 and rf = (t* ±irj)/2 
denote the raising/lowering operators for the two spin 
species] . 

The spin current j analogous to the electron current 
(I13p is derived from the continuity equation for the local 
spin operators S^ = (j%/2: dtS^+'V(j (T )k = 0, which can 
be rewritten as (j a )k+i ~ (ja)k = ^[<7 z k ,H s ] (analogous 
equations can be written for the r species in Eq. (jT4j) ) . 
We obtain 



J 

jo 


= 3a +Jr, 




jr 




1 T k T k+l>- 



In the following, we choose a symmetric driving: Tl = 
T R = T and f LtR = |(1 T /), so that / = f R - f L G [0, 1] 
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(/l < /r, < f L ,fn < 1) is the parameter controlling 
the driving strength. Small / implies that the system 
is weakly driven by the external baths, and behaves as 
in the linear response regime. In the opposite limiting 
case / = 1, the left (right) bath only induces up-down 
(down-up) spin flips for both spin species. 

Using the method of quantum trajectories we evalu- 
ated the stationary spin currents (j T ) for the two 
species of spins. Due to the mapping between the elec- 
trons described by the Hamiltonian in Eq. (p~2|) and the 
spins obeying Eq. (|14p . this spin current exactly equals 
the electronic current in the Hubbard model. In particu- 
lar, (j a ) ((J T )) is the current flow of electrons with spins 
pointing up (down), that is the crucial physical quantity 
in charge transport. 

Perhaps the most interesting result we found in the 
current behavior as a function of the driving is the emer- 
gence of a NDC phenomenon for sufficiently strong driv- 
ings, as shown in Fig.O It happens that, while for small 
/ values the current increases, there exists a value /* at 
which (j) exhibits a maximum and then, further increas- 
ing /, it decreases. 

One can now question whether or not this non- 
monotonic behavior is stable when varying the Hamil- 
tonian parameters t and U . In all simulations reported 
here we fixed energy units by setting t — 1. In the lim- 
iting case where U = 0, the fermions in the Hubbard 
model are non-interacting, therefore a linear regime in 
which the current is always proportional to the driving 
strength is expected. In view of these considerations, it 
is tempting to assume the existence of a critical value 
U* in the Hamiltonian parameters space separating the 
linear and NDC behaviors of the current. As a matter of 
fact, within numerically accessible system sizes, we ob- 
served NDC only for U > U* , where U* « 2. This can 
be seen from Fig. [3J where we plot the maximal current 
drop, measured by (j) / = /» — (j) f—i, as a function of the 
on-site interaction strength U. From the inset it is clear 
that, while for U <C 2 the current is proportional to the 
driving, for U 3> 2 a bell-shaped behavior emerges. Of 
course, on the basis of the data presented in Fig. [3J one 
cannot exclude that U* drops with N. In this scenario, 
in the thermodynamic limit NDC would be observed for 
any U > 0; nonetheless, we point out that a mean-field 
qualitative argument given at the end of Sec. IV Dl sup- 
ports the existence of NDC for U > 2, thus agreeing with 
our findings in Fig. [3] In any case, a significant result of 
our numerical simulations is the emergence of NDC in a 
physically relevant transport model such as the Hubbard 
model, at small system sizes N > 4. 



IV. SPINLESS FERMION MODEL 

The investigation of the far-from-equilibrium proper- 
ties of the Hubbard model is numerically demanding and 
an analytical treatment appears difficult. Therefore, in 
what follows we will focus on a simplified model, usually 
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FIG. 3: (Color online). Maximum current minus current at 
maximum driving strength, (j)/* — (j)i as a function of the 
on-site repulsion U. In the inset we show the spin current as 
a function of /, for a fixed size N = 6 and different values 
of U: from top to bottom U — (circles), 0.5 (squares), 1 
(diamonds), 1.5 (triangles up), 2 (triangles left), 2.5 (triangles 
down), 3 (triangles right), 5 (stars). Data are for F = 1, and 
the simulation time (QT approach) is T = 10 5 . The current is 
plotted only for the er-spin species; differences with the r-spin 
species are negligible in the scales of the figure. 



referred to as the t — V model, where the spin degree of 
freedom is neglected. This model is numerically much 
more convenient. Moreover, it will help us in gaining a 
deeper understanding of the peculiarities of charge trans- 
port discussed in Sec. IIIII for the Hubbard model. 

The t — V model considers spinless fermions instead of 
spin-1/2 particles. Its Hamiltonian reads as follows: 



c ]+i c j) 



(16) 



Similarly to the Hubbard model (|12p , the operators cj , Cj 
create/annihilate a spinless fermion at site j = 1, . . . , N 
(therefore they satisfy canonical anticommutation rules) , 
while rij — CjCj is the corresponding number operator. 
The system parameters t and V describe, respectively, 
the nearest neighbor hopping strength and the fermionic 
repulsion between contiguous sites. 

In direct analogy with what has been discussed for the 
Hubbard model, we take open boundary conditions and 
couple both ends of the chain to some external baths that 
inject and extract fermions. Now the number of Lindblad 
operators is halved, since we removed the spin degree of 
freedom: 



L3 = a/TrTr ct 



(17) 



-JV > 



L 2 = y/T L (l-f L )c 1 , 
La = y/r R (l - f R ) c N ; 

the parameters T^, T^?, and fu play roles analogous 
to those of the corresponding parameters introduced for 
spinful fermions in Eqs. ([5]), (|6"]). 
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The Hamiltonian in Eq. ([16]) can be mapped into an 
XXZ Hciscnbcrg spin chain plus some spurious contribu- 
tions, consisting in external transverse magnetic fields, 
which arc irrelevant in the isolated (Hamiltonian) case 
and do not qualitatively modify the spin current behav- 
ior (see Appendix IB). Therefore we shall neglect these 
spurious terms and concentrate on the XXZ spin-1/2 sys- 
tem whose autonomous Hamiltonian is given by: 



n, 



N-l 

E 

i=i 



[u 



(18) 



where (a — x,y,z) are the Pauli matrices of the j- 
th spin, and A = J z /J x denotes the xz anisotropy; N 
is the total number of spins. In the t — V model lan- 
guage of Eq. (|16p . the couplings in Eq. (fl8|) are given 
by J x = —t/2 and J z — V/A, so that the anisotropy 
A = —V/2t. Strictly speaking, since fermionic inter- 
action between contiguous sites is repulsive (V > 0), 
this would correspond to antifcrromagnetic transverse 
couplings J z > 0. Nonetheless, when considering the 
XXZ spin model (fT5|) . one is not a priori forced by this 
constraint and can also analyze the ferromagnetic case 
J z < 0. Hereafter we set J x = 1 as the system's energy 
scale. 

In the spin-1/2 picture, the Lindblad operators (JT7J) 
are mapped into operators flipping the border spins. In- 
deed as explained in Appendix [Bj apart from uninfluent 
phase factors, we have c\ — > af and c, — > a~ . The 
Fermi function Jl.r is such that — Is [—1,1] is 

the corresponding bath's magnetization per spin in di- 
mensionless units. As we did for the Hubbard model, we 
choose (with the exception of Sec. IVB4[) a symmetric 
driving: T L = T R = V and f LM = ±(1 qp /), so that 
f = f R -he [0, 1] (f L < f R , 0<f L J R < 1) is a single 
parameter controlling the driving strength. When / is 
small we are in the linear response regime, while in the 
limiting case / = 1 (corresponding to Jl = 0, fu = 1) 
the left (right) bath only induces up-down (down-up) spin 
flips. The spin current is computed as in Eq. (|15p. but 
without the contribution of the r species: 



(19) 



Quantitative numerical and semi-analytical analysis of 
the model (fT5)) are easier than in the model p4[) . In par- 
ticular, the local Hilbert space is halved: for a fixed num- 
ber of sites N, the size Af = 2 N of a generic state vector 
describing the system is decreased by a square root factor 
with respect to the size (2 N ) 2 of a state vector for spin 
ladder (|14[) of length N. It is therefore clear that, with- 
out truncating the Hilbert space, using the Monte-Carlo 
wave function method one is able to simulate chains of 
twice the length of a ladder with the same computational 
cost. We wish to note that in the spin ladder systems we 
have so far used only the QT method to perform nu- 
merical simulations. Also the MPO approach could in 
principle be used, by simply joining two sites from the 



opposite spin chains into a single site with a local dimen- 
sion 16 (for a single chain it is 4). The complexity of 
time evolution increases by a factor 4 3 at a fixed matrix 
dimension D (see Appendix lAl for definitions), because 
of singular value decompositions of 16D x 16D matrices, 
instead of AD x AD. 



V. SPIN TRANSPORT PROPERTIES IN A 
HEISENBERG CHAIN 

In this section, we study the far-from-equilibrium 
transport properties of the XXZ Heisenberg spin 
chain (fT5)l , with the two edge spins coupled to exter- 
nal baths, as described in the previous section. We first 
present, in Sees. IV AIIV CI the results of our numerical 
simulations, focusing on the NDC phenomenon and on 
the appearance of a long-range spin ordering into ferro- 
magnetic domains. Then, in Sec. IV Dl we qualitatively 
explain our results in terms of one-magnon localization. 
Since our findings are suggestive of a phase transition 
with the emergence of long-range order, we have also 
searched for numerical evidence of such transition by an- 
alyzing the spin-spin correlation function. Our data dis- 
played in Appendix [U] even if not conclusive, show a 
dramatic slowing down of the correlation decay in the 
NDC regime, even though the accessible system sizes are 
too small for a quantitative analysis of a possible phase 
transition. Finally, in Sec lVEI we rephrase the results in 
terms of the fermionic current. 

Note that the spin current in the XXZ model with anti- 
ferromagnetic coupling is the same as the charge current 
in the t— V model, therefore all the results we discuss here 
for J z > also apply to the case of fermionic transport in 
Eq. (Til))) . For the spin transport we also considered cases 
where J z < 0, thus corresponding to a rather unphysi- 
cal attractive fermionic interaction in the t — V model; 
quite surprisingly, we found that data obtained for the 
XXZ chain are insensitive to the sign of the transverse 
coupling J z . 

As far as we know, transport properties of the 
autonomous model described by Eq. (|18p have been 
extensively analyzed only within the linear response 
even if a fully comprehensive under- 
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regime- 
standing is still lacking. In particular, it has been found 
that the low-temperature and the high-temperature ther- 
modynamic transport properties are essentially deter- 
mined by the xz anisotropy Ai£. In the zero magne- 
tization sector, the XXZ model is an ideal conductor for 
|A| < 1, while numerical dataii suggest that the sys- 
tem is a normal (diffusive) spin conductor for |A| > 1. 
The normal conduction in the |A| > 1 regime has been 
recently confirmed for systems of much larger size^, 
N ~ 100, see also Ref.[H. The above two distinct behav- 
iors may be associated to two different system phases at 
zero temperature: for — 1 < A < 1 the system is gapless, 
while for A < —1 (A > 1) it is ferromagnetically (antifer- 
romagnetically) ordered, and the ground state exhibits a 



finite gap with the first excited state. 

We now investigate the far-from-equilibrium proper- 
ties of the XXZ Hcisenberg spin chain, beyond the linear 
response regime. 



A. Gapless regime 

We start by considering the gapless phase, where the 
linear response theory predicts ballistic transport. We 
found that in the regime |A| < 1 the current is always 
proportional to the driving and independent of the chain 
length TV; this holds for any value of the driving strength 
/ G [0,1]. Remarkably, there are no appreciable quan- 
titative differences between the current evaluated with 
a ferromagnetic (J z < 0) or antiferromagnetic (J z > 0) 
coupling. Fig. [4] displays the spin current as a function of 
the driving strength /, for A = 0.5. The corresponding 
spin magnetization profiles for two distinct values of / 
are displayed in the insets as a function of the site in- 
dex. In both cases we can see a nearly flat profile, that is 
typical of systems with ballistic spin propagation^ 3 -; most 
importantly, we notice that the stationary spin magne- 
tization at the borders is very different from the bath 
magnetizations (o z L r ) = =F./\ 



N(j) 





-*N= 


4 




-•N = 


6 




N = 


8 




-* N= 


12 


T" 


*N = 


16 




-<N = 


40 




FIG. 5: (Color online). Spin current as a function of the 
driving strength /, in a chain with anisotropy A = 2. The 
system-bath coupling is set equal to T = 4. Symbols for 
TV < 16 are obtained using QT after a simulation time T = 
2.5 x 10 3 , while curves display data for a longer integration 
time T = 7.5 x 10 4 . Data for TV = 40 are obtained by MPO, 
with integration times up to T = 4000. Long-time data for 
TV < 16 are quoted from Ref. OS- 




FIG. 4: (Color online). Spin current as a function of the 
driving strength for A = 0.5. The system-bath coupling is set 
equal to V = 1. The insets show spin magnetization profiles 
versus the scaled spin index coordinate, for two values of / = 
0.5, 1. Data are obtained from the QT approach. 



B. Gapped regime 

The transport properties are much more interesting in 
the gapped phase where, in the linear response regime 
/ <C 1, numerical data suggest normal diffusive trans- 
por t 17 ! 33 . For | A| > 1 the spin current is no longer mono- 



tonic with / and exhibits a typical bell-shaped behavior, 
as we observed for the Hubbard model (see Fig. \2§ ■ 

In Fig.[5]we plot the spin current as a function of /. For 
small / the system behaves as a normal Ohmic conduc- 
tor, as expected from linear response result a 17 i 33 : namely 
the current increases like (J) oc / /TV. After a given value 
/* of the driving at which the current reaches its maxi- 
mum, it starts decreasing with / until it is strongly sup- 
pressed at / = 1. Details on the scaling of (J) with the 
system size at / = 1 are given in Sec. lVBTl Here we just 
point out that, interestingly, NDC is already visible af- 
ter short integration times: as a matter of fact, with the 
QT approach (see data for TV < 16), the characteristic 
bell-shaped behaviors can be obtained after a simulation 
time T ~ 2.5 x 10 3 (symbols) that is much shorter than 
the one required to reach the stationary values, T ~ 10 s 
(curves); furthermore, NDC features are present even af- 
ter very short times T~ 5 x 10 2 . 

The data shown in Fig. \5\ are suggestive of a transi- 
tion from a normal spin conductor phase at small /, to 
an insulator phase at large /. On the other hand, for 
the achievable system sizes, the value f*(N) where the 
current reaches a maximum drops with TV. Therefore, 
in principle we cannot exclude an alternative scenario 
where at the thermodynamic limit the system becomes 
an insulator at any driving strength /. In any case, as 
far as the bias / is increased, a substantial modification 
of spin transport properties becomes apparent. These 
results are suggestive of a far-from-equilibrium quantum 
phase transition. In the light of a recent paper of one 
of us^ 2 -, we have analyzed the spin-spin correlation func- 
tions, in order to see the possible emergence of a phase 
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FIG. 6: (Color online). Spin magnetization profiles versus 
scaled spin index coordinate for the same parameters as in 
Fig. [5] and driving strengths / = 0.3 (empty symbols) and 
f — 1 (filled symbols). Data are obtained from the QT ap- 
proach. — Figure acknowledged from Rcf. 




FIG. 7: (Color online). Spin current at a fixed driving 
strength / = 1 for the Heisenberg spin model with A = 2. 
The system-bath coupling is set equal to T = 4; the simula- 
tion time (QT approach) is T — 7.5 x 10 4 . The thick curve 
displays an exponential fit of numerical data (circles). 



transition that should be characterized by the emergence 
at strong driving strength / of a long-range correlation 
order. Numerical data, even if not conclusive, are shown 
in Appendix [Cl and are in support of such a behavior. 



1. Behavior at f = 1 

As hinted in Ref. [32], in order to understand the physi- 
cal mechanism lying at the basis of NDC, we have to an- 
alyze the stationary spin magnetization profiles. These 
are shown in Fig. [5] Note that, in contrast with the fast- 
time raising up of the NDC phenomenon, a much longer 
integration time is required in order to reach a good con- 
vergence for the spin magnetizations, due to the equili- 
bration time scales that, at / = 1, grow exponentially 
with the distance of the spin from the chain border—. 

As shown in Fig. [SJ magnetization profiles in the lin- 
ear response regime / < 1 exhibit a constant linear gra- 
dient, where the magnetizations of the two edge spins 
are close to the bath magnetizations =p/, as it is ex- 
pected for normal Ohmic conductors. In the limiting 
case / = 1 a peculiar stationary state characterized by 
two ferromagnetic domains that are oppositely polarized 
appears. Moreover their relative width increases with 
the system size. These domains are eventually responsi- 
ble for strongly suppressing the spin current, since they 
inhibit spin flips. Strictly speaking, evidence of the for- 
mation of such domains is also visible for smaller, though 
strong drivings, where the NDC effect is established (see 
Sec. IV B 2|) . We will explain later in Sec. lVDl the physical 
mechanism leading to the formation of such domains in 
the gapped phase. Here we stress that, as for the gapless 
regime, we found no quantitative differences between a 



ferromagnetic or antiferromagnetic coupling. This may 
be somewhat counterintuitive, since in this last case fer- 
romagnetic domains correspond to a highly excited state 
for the autonomous system; the mechanism therefore has 
its roots in the genuine far-from-equilibrium dynamics. 

We already observed that the current for small driv- 
ing strengths behaves as for a normal Ohmic conductor: 
(?) f/N. On the other hand, Fig. [7| shows that at 
maximum bias / = 1 the current drops to zero exponen- 
tially with N, that is, the system is an insulator. Un- 
fortunately we could not achieve system sizes larger than 
N = 12, since there the current is very small, therefore 
it would require a huge number of time steps in order to 
get reliable results, thus making simulations unfeasible 
(already for N — 12 the simulation time T — 7.5 x 10 4 
is not long enough: the last point in the figure has been 
obtained by observing only one spin flip during the whole 
QT simulation). 

2. Spin blockade 

At maximum driving the current is strongly sup- 
pressed, due to the fact that ferromagnetic domains of 
macroscopic length ~ N/2 are formed. Nonetheless, a 
strong inhibition of the spin current can be also achieved 
by only creating a much smaller ferromagnetic region 
close to each bath. Indeed, signatures of this spin block- 
ade mechanism are already seen at / ~ 0.9 — 0.95, where 
asymptotically only a couple of outer spins reach mag- 
netization values close to ±1, but still the current is far 
below its maximum value (j) /*. 

Some spin magnetization profiles for strong drivings 
are explicitly shown in Fig. [5J The dotted-dashed green 
curve corresponds to a maximum driving; there macro- 
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FIG. 8: (Color online). Spin magnetization profiles versus 
scaled spin index coordinate in the insulating case A = 2 
with a system-bath coupling T — 4, for strong drivings. In the 
main panel we fix a system size N = 12 and choose different 
values of /, while in the inset we focus on the case / = 0.9 
and vary the size according to the legend. Data are obtained 
from the QT approach. 
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FIG. 9: (Color online). Spin magnetization profiles versus 
scaled spin index coordinate for a chain with N — 8 spins 
at maximum driving (/ = 1), for different values of A > 1; 
the system-bath coupling is set equal to F = 1. In the inset 
we plot an estimate of the thickness of the interface region 
£ as a function of hi(A); the straight line shows the fit £ tx 
(In A) -1 ' 625 . Data are obtained from the QT approach. 



scopic ferromagnetic domains are clearly visible. The 
other ones are for / slightly less than one, but it is still 
possible to see that a couple of spins close to the bor- 
ders are nearly perfectly down/up polarized. In the inset 
we fix / = 0.9 and vary the system size; we notice that, 
when increasing N, the number of spins involved in the 
spin blockade also increases, in accordance with the re- 
sults previously shown for / = 1. 



3. Thickness of the interface region 

An interesting point that can be addressed is the anal- 
ysis, at maximum driving strength / = 1, of the charac- 
teristic thickness £ of the interface region, located around 
the chain center and dividing the two ferromagnetic re- 
gions. According to our data shown in Fig. [9j this size 
depends on the system anisotropy A. On the other hand, 
we checked that dependence on N is negligible. In order 
to give a rough estimate of £, we fixed a threshold (crj) 
(horizontal dot-dashed line in the figure, corresponding 
to (tr|) = 0.6) and then evaluated the distance £ of the 
point in each spin magnetization profile reaching that 
value from the limiting case (A — > oo) in which the chain 
is exactly split into two ferromagnetic domains (due to 
the fact that our problem is on a lattice of finite length, 
the A — > oo magnetization for N/2 < j < N/2 + 1 is 
estimated after joining the two perfectly ferromagnetic 
domains with a skew dashed line). In the inset we study 
the dependence of the distance £ on A; the straight line 
shows the fit £ oc (In A) -1 - 625 . We point out that the 
semi- analytical argument of the one-magnon localization 



length that will be discussed in Sec. IV Dl predicts that the 
size £ of the interface region is of the order of the one- 
magnon localization length, namely logarithmic in A and 
TV- independent . 

4- Stability of the ferromagnetic domains 

Until now we have only considered situations in 
which system-reservoir couplings Tl.r and the driving 
strengths fiL,R are symmetric. At this step one may won- 
der if the position of the domain wall between the two 
ferromagnetic regions in the spin magnetization profiles 
is stable against the breaking of such symmetry. Any im- 
balance, though small, may in principle cause a shift of 
the interface region towards one of the boundaries of the 
chain. This would raise some doubts about the stability 
of the previously depicted scenario, making our discus- 
sion relevant only for fine-tuned values of the Lindblad 
parameters. Below we show that this is not the case. 

Imbalances of the couplings Tl^r have quite tiny ef- 
fects on the steady state at maximum driving, as one 
can see from the upper panel of Fig. [TO] We put there 
a strong asymmetry, by setting Tl = 1, =0.1, and 
thus admitting a coupling to the right bath that is one 
order of magnitude smaller than the one to the left bath. 
Differences with respect to the symmetric case are ap- 
parent for finite integration times T, where the domain 
wall is clearly shifted to the right. On the other hand, 
as far as T is increased, profiles become more symmetric. 
It is not a priori clear whether the steady state is per- 
fectly symmetric, as those in Fig. [6] from our data we 
cannot rule out possible deviations in the position of the 
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FIG. 10: (Color online). Spin magnetization profiles versus 
spin index coordinate for a chain with TV = 10 spins and 
anisotropy A = 2, at different integration times T. Upper 
panel: symmetric drivings such that / = 1 and system bath 
couplings r_L = 1, r_R = 0.1. Lower panel: the system-bath 
coupling for the left and the right baths is kept fixed and equal 
to r = 1, while the drivings are chosen asymmetrically as 
Sl = 0, fit — 0.99. Data are obtained from the QT approach. 



domain wall which are logarithmic in the coupling im- 
balance; this would be hardly detectable from a merely 
numerical analysis. 

Stronger modifications are induced by imbalances on 
the driving strengths (J-l.r- Indeed in that case even a 
small asymmetry would cause a weaker spin blockade on 
one side, as discussed in Sec. IV B 2l (see the lower panel of 
Fig. [TU1 in which we put Jl = 0, Jr = 0.99); as a conse- 
quence, a deformation and a broadening of the interface 
region are also established. Nonetheless, we point out 
that these modifications appear to be continuous in the 
degree of imbalance, and thus in principle controllable. 



C. XXX Heisenberg isotropic model 

The isotropic XXX Heisenberg chain (that is, |A| = 1) 
corresponds to a limiting, nonetheless interesting situa- 
tion, since molecular compounds that are used to inves- 
tigate one dimensional spin- 1/2 transport properties are 
often very well described by isotropic antiferromagnetic 
Heisenberg exchange coupling o 39 i 54 . 

In Fig. [Tl] we show some numerical data concerning 
the behavior of the spin current with respect to the driv- 
ing field: NDC is visible only for sufficiently long chains 
(TV > 8). A qualitative understanding of this result 
comes from an analysis of the spin magnetization pro- 
files at / = 1, that are plotted in the inset. As a matter 
of fact, we can recognize a situation that is similar to the 
one already observed at |A| > 1 and / < 1 (see Fig. [SJ: 




FIG. 11: (Color online). Spin current as a function of the 
driving strength for the isotropic case A = 1. The system- 
bath coupling is set equal to T = 2. The simulation time 
(QT approach) is T = 5 x 10 4 . In the inset we show the spin 
magnetization profiles versus scaled spin index coordinate, at 
maximum driving / = 1. Different curves stand for various 
system sizes. 



by increasing TV, a partial spin blockade of the outermost 
spins is established. This progressively inhibits the cur- 
rent flow along the chain, thus setting up the NDC mech- 
anism. Notice also that, at small TV values (TV = 4, 6), 
the spin blockade is very weak; this prevents systems of 
very small size from exhibiting the NDC phenomenon, 
even though a nonlinear dependence of the spin current 
on the driving strength can already be seen. 

We conjecture that also for the isotropic XXX chain 
the NDC phenomenon is stable at the thermodynamic 
limit. Indeed, our simulations suggest the following pic- 
ture: on one hand, at small / the current decreases as 
(j) oc TV~ a with a « 0.4 (data for / = 0.2 are shown 
in Fig. rjJl left; we checked that this is consistent for all 
/ 0.3). On the other hand, at / = 1 it drops to zero 
faster than linearly with TV (probably exponentially, as 
suggested by Fig. [T21 right), thus indicating a relative 
current drop 1 — (j)/ = i/(j)/ = /« that increases with TV 
towards the unit value; here /* denotes the driving at 
which the current is maximal. In both panels a log-log 
scale has been used, so to make visible the distinction be- 
tween a power-law scaling at small / and an exponential 
behavior at / = 1. We also plotted the 1/TV behavior 
(dashed lines) expected for normal Ohmic conductors, 
such to show that for / <C 1 the decay is slower than 
that, and for / = 1 it is faster). 



D. One-magnon localization 

The formation of ferromagnetic domains and the nega- 
tive differential conductivity phenomenon in the gapped 
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FIG. 12: (Color online). Spin current for A = 1 at a fixed 
driving strength, / = 0.2 on the left, while / = 1 on the right 
(see Fig. Ill|l . Full curves are fits of numerical data (obtained 
from QT): while an inverse power-law fit works well at / <C 1, 
for maximum driving an exponential fit seems adequate. The 
dashed lines indicate a behavior (j)/ ~ 1/N and are plotted 
as guidelines. 



regime | A| > 1 can be qualitatively explained in terms of 
localization of one-magnon excitations created at the bor- 
ders of a ferromagnetic domain. We should immediately 
point out that the following argument is independent of 
the ferromagnetic of antiferromagnetic coupling. This 
reflects into the fact that the steady state with ferromag- 
netic domains is completely driven by dynamical effects, 
and not by the system ground state, being actually an 
antiferromagnet for J z > 1. 

Given a ferromagnetic state |0) = \[[ •••!), one- 
magnon excitations have the general form X^fcLi a k 
where | k) = |0) describes the state with the fc-th spin 
flipped. If the autonomous XXZ chain has open bound- 
ary conditions, there is an energy gap 2\J Z \ between the 
states |1) and \N) (spin- flip excitations at the bound- 
aries) and the states |2), |3),..., \N — 1). Indeed, we have 

(0\H S \0) = (N-1)J Z , 

(l\H,\l) = {N\H,\N) = (N-3)J m , (20) 
<2|W fl |2) = ... = (N - 1\H S \N -1) = (N- 5)J„ 

where TC S is the XXZ Hamiltonian (fP5)) . Only nearest- 
neighbor spin-flipped states are coupled and the coupling 
strength is 2\J X \: 

(k\H s \k + l) =2J X , k=l,...,N-l. (21) 

As shown in Appendix [Dj the autonomous model (12"0)) - 
(f2"Tj) is exactly solvable in the limit of large N and spin- 
flip excitations created at the borders of the chain remain 
exponentially localized when | J z \/\ J x \ — |A| > 1, over a 
localization length I ~ 1/ In |A|. 

We now consider the coupling to external baths. First, 
it is instructive to discuss the case in which the system 
is coupled to a single, fully polarized reservoir, = 0. 
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FIG. 13: (Color online). Spin magnetization profiles at A = 
2, r = 4 for a maximal driving strength: / = 1. The various 
panels correspond to different times T (in QT simulation) at 
which the profiles are plotted. Circles are for a single bath 
coupled to the first spin (j — 1), while squares are for two 
baths at the ends of the chain. 



Regardless of the anisotropy A, the stationary state is 
pure and ferromagnetic, namely \H ■ ■ ■ J.) (J, J, ■ • • J.|, since 
the Hamiltonian (|18p conserves the overall magnetization 
while at the left boundary of the chain only the lower- 
ing operator L2 oc crj~ acts (note that the convergence 
to the stationary ferromagnetic state can be rigorously 
proven following Ref. |25T ). As shown in Fig. [13] (see 
circles), the time scale required for the convergence of 
the j-th spin to the equilibrium state |J.) scales exponen- 
tially with j. Consider now an intermediate state with 
m spins down. To enlarge the ferromagnetic domain, 
one-magnon excitations should be propagated, through 
<7j(Tj +1 and crj(rj' +1 exchange couplings of the Hamilto- 
nian (fT8|) . across the ferromagnetic domain to the left 
chain boundary. Suppose, for instance, that we have the 
leftmost m spins down and the (m + l)-th spin up, and 
that this excitation propagates to the left bath; then the 
bath can flip this spin down, thus ending up with a fer- 
romagnetic domain with m + 1 leftmost spins down. The 
crucial point is that the one-magnon propagation is ex- 
ponentially localized at |A| > 1. Hence, exponentially 
long time scales cx exp(j/£) ~ exp(jln|A|) are required 
to polarize the j-th spin. 

Consider now two baths at / = 1. Due to exponen- 
tial one-magnon localization, essentially only the nearest 
bath is felt by the spins in the chain, with the exception 
of an interface region between the two ferromagnetic do- 
mains, whose length is of the order of the one-magnon 
localization length. As shown in Fig. [13] (squares), apart 
from the interface region, the spin magnetization profiles 
for the spins closer to the left than to the right bath 
essentially evolve as in the single-bath case. 



12 



These ferromagnetic domains are responsible for 
strongly inhibiting spin flips, and therefore for suppress- 
ing the spin current at / = 1. Since at small / the 
current grows linearly, we can conclude that, due to the 
continuity of (j)f, a region of negative differential con- 
ductivity exists. Finally, we note that, in agreement with 
our numerical data, the one-magnon argument does not 
distinguish between ferromagnetic ( J z < 0) and an anti- 
ferromagnetic (J z > 0) spin couplings. 

The argument developed in this section can be ex- 
tended to the spin ladder model (|T4]) obtained by apply- 
ing a Jordan Wigner Transformation (JWT) to the Hub- 
bard Hamiltonian. Assume that we have a ferromagnetic 
domain of the to spins of both species (a and r) clos- 
est to, say, the left bath. Now consider a spin flip for 
the (to + l)-th spin, for instance of the a species and as- 
sume that the spins of the r species can be treated within 
the mean-field approximation. The only energy term not 
constant for the autonomous evolution (fT4f restricted to 
the one-magnon sector (a species) is ^ . a j T j- Within 



mean- field approximation, we substitute tJ — > (tJ), and 

= 0. There is 



= 1 for 



1, m, while (t. 



m+l/ 



an energy gap ^ between states with the spin flipped 
belonging to sites from 1 to to and the state with the 
(to + l)-th spin flipped. The hopping strength is \t\. 
Therefore magnon localization as in the XXZ model takes 
place for \U\ > 2\t\. Such prediction is compatible with 
the numerical results shown in Fig. [3] where NDC for 
the Hubbard model is observed only for U > U* , with 
U* w 2. However, the argument is of a mean field nature 
and therefore has to be considered weaker than the one 
developed for the XXZ chain. 



E. Discussion in terms of the charge current in the 
fermionic model 



At this stage, it is useful to summarize the main re- 
sults obtained for the spin chain in terms of the orig- 
inal electronic t — V model. Since |(1 + (<r^)) corre- 
sponds to the electronic charge density at site k, the fer- 
romagnetic domains observed at / = 1 corresponds, in 
the fermionic picture, to a phase separation, with all the 
electrons frozen in the right half of the lattice, close to 
the emitter electrode. Therefore, charge transport is in- 
hibited, provided that V/t > 2. Note that this charge 
clustering takes place in spite of the repulsive nature of 
electron-electron interactions. The magnon localization 
argument can be straightforwardly reformulated in terms 
of the fermionic model. In particular, a one-magnon ex- 
citation becomes a single electron (hole) propagating on 
an empty (filled) lattice. The single-bath case can then 
be interpreted in terms of depletion (filling) of the lattice 
by means of a single lead, playing the role of a charge 
collector (emitter). This process requires exponentially 
long time scales at V/t > 2. Finally, we point out that 
the NDC regime is observed for strongly interacting sys- 
tems (V/t > 2) and in the far-from-equilibrium regime, 



corresponding to large bias voltages eV ^> fcgT, so that 
the Fermi functions for the collector and the emitter elec- 
trodes, evaluated at the energy differences E\ and En, 
satisfy /l ~ and /jj k 1, respectively. 



VI. NONINTEGRABLE MODEL: STAGGERED 
MAGNETIC FIELD 

The high-temperature transport properties in one- 
dimensional quantum many-body systems are strongly 
affected by the presence of conservation law a 20 ' 21 i 55 ' 56 . 
In particular, the existence of local conserved quantities 
Q n , n = 1, 2, . . ., typically leads to an ideally conducting 
- ballistic behavior, at all the temperatures. This is a con- 
sequence of the inequality due to Mazur— which bounds 
the time averaged current-current autocorrelation func- 
tion as hm t ^ 00 (l/t) f*dt'(J(t')J(O)) > Y, n \( J Qn)p? 
where (X)p = tr [exp(-/3H)X]/ tr [exp(-/3H)] and Q n 
are chosen and normalized such that (Q n Q m ) — 5 nm . 
As one can see, this argument essentially depends on, 
first, the existence of nontrivial conserved quantities (as 
is typically the case only for completely integrable sys- 
tems), and, second, on the overlaps (Q n J)p between the 
conservation laws Q n and the transporting current J in 
question. For example, for the XXZ model all the con- 
servation laws have zero overlaps with the spin current, 
(J s Qn) = 0j which allows in the gapped regime |A| > 1 
for the normal diffusive spin transport as discussed ear- 
lier, whereas, for example, the heat transport is ballistic 
as the energy current is just one of the conserved quanti- 
ties Je = However, even though the XXZ Heisen- 
berg model is completely integrable, we have seen that 
the spin transport properties in the far-from-equilibrium 
regime are very different from the linear response regime 
behavior. Therefore, we might expect that the presence 
of NDC is not related to integrability of the Heisenberg 
model. 

In order to check the stability of the nonlinear trans- 
port highlighted in the previous sections with respect to 
breaking system integrability, we considered a slightly 
modified spin chain model, in which a staggered mag- 
netic field along the z direction is added to the Heisen- 
berg Hamiltonian of Eq. (|18|) . Namely, we studied the 
following autonomous model: 



Hs 



[w+i- 



j u 3 + l 



N 



■lVo-2 



(22) 

Interestingly, the model (f2"2"]l exhibits a transition from 
integrability to quantum chaos when increasing the field 
strength B. This can be detected in the change of the 
spectral statistics of the system 5 ^' 59 . In particular, in 
Fig. [TJ] we plot the integrated level spacing distribution 
I(S) [I(S) is the probability that a randomly chosen level 
spacing - normalized to the mean level spacing - is less 
than S]. It is shown that, for a given set of parame- 
ters and B ^ 0, the level statistics follows the universal 
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FIG. 14: (Color online). Integrated level spacing distribution 
I(S) for the Hamiltonian in Eq. J22|) with A = 1.3, B = 0.3, 
N = 16; data (full black curve) correspond to the zero magne- 
tization sector. In order to apply the random matrix theory, 
the system Hamiltonian has to be diagonalized in a subspace 
in which no symmetries (except the time-reversal) are left. 
For this purpose, in addition to fixing the total magnetiza- 
tion, we applied the staggered magnetic field on all spins but 
the first one (i.e., we supposed that the last sum in Eq. (|22[) 
runs from 2 to N). This breaks the spatial reflection sym- 
metry j — ► N — j, without affecting the transport properties 
under investigation. Red dashed (blue dotted) curve indicates 
Wigner-Dyson (Poissonian) statistics, typical for chaotic (in- 
tegrable) systems. In the inset (data from QT) we plot the 
spin current as a function of the driving strength, for different 
chain lengths; we fixed a system-bath coupling V = 1. 



predictions of the random matrix theory (Wigner-Dyson 
statistics in presence of time-reversal symmetry) , as typ- 
ical for chaotic (strongly non-integrable) systems. In the 
inset we analyzed the corresponding spin current as a 
function of the driving strength: we found a qualitatively 
analogous behavior as in the integrable case, with a NDC 
regime still clearly visible. We also checked the presence 
of a normal Ohmic conduction for small drivings: from 
Fig. H~5l one can see that, for / = 0.1, the spin current 
scales as (j) oc 1/N, according to Ohm's law of diffu- 
sive transport (see also Ref. l33l). therefore spin transport 
is normally diffusive, as expected for a chaotic system. 
In analogy with the integrable case, for larger gradients 
where the current behaves highly nonlinearly, the spin 
current decays faster than linearly (at / = 0.5 we found 
a decay (j) ~ AS/N 1A , thus indicating an insulating 
behavior in the thermodynamic limit). Also magneti- 
zation profiles in that case are not linear anymore, and 
the emergence of a weak spin blockade is already visible 
at / — 0.5, for sufficiently large sizes (see the inset of 
Fig. I15p , similarly to what has been already discussed in 
Sec. ESI 

On the other hand, we considered a situation in which 
the system does not exhibit NDC (A = 0.5, B = 0.5), 



FIG. 15: (Color online). Dependence of the scaled spin cur- 
rent on the chain length in the staggered Heisenberg model, 
with A = 1.3, B — 0.3, and F = 1 (AS is the magnetiza- 
tion difference across the chain), at driving strength / = 0.1 
(circles) and / = 0.5 (squares). The two dashed lines re- 
spectively denote the behaviors ~ 1.3/N ' , (for circles) and 
~ 1.85/-/V 1 ' 4 (for squares), whereas dotted line indicates 1/N 
behavior found in the linear response regime. In the inset we 
show magnetization profiles for / = 0.5 and different system 
sizes. Data are obtained from MPO. 



but still it is quantum chaotic with respect to energy level 
statistics. In that case we found the usually predicted 
behavior for non integrable systems^: the current is al- 
ways proportional to the driving for any value of /, and a 
normal Ohmic regime for both small and strong drivings 
emerges, as shown in Fig. 1161 This comes in sharp con- 
trast to the integrable case, where at A < 1 ballistic spin 
transport takes place. Consistently with normal metal- 
lic conductors, the spin magnetization profiles exhibit a 
linear gradient (see the inset of Fig. [T6|) . 

From the results of this section we can conclude that 
the NDC phenomenon is governed by the anisotropy pa- 
rameter A and is insensitive to the transition from in- 
tegrability to quantum chaos. The appearance of NDC 
in the quantum chaos regime is particularly significant 
since in this case normal Ohmic conduction is expected 
in the linear response regime (small driving /). Due to 
the insulating behavior observed for |A| > 1 at large 
driving strengths /, the system should undergo a metal- 
to-insulator transition when increasing /. 



VII. CONCLUSIONS 

In summary, we have performed extensive numerical 
simulations showing that the one dimensional Hubbard 
model, as well as the corresponding spinless fermion 
model, exhibits a regime of negative differential conduc- 
tivity, in which the particle current decreases as one in- 
creases the driving field. The same phenomenon, trans- 
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FIG. 16: (Color online). Dependence of the scaled spin cur- 
rent on the chain length for the staggered Heisenberg model, 
with A = 0.5, B — 0.5, and V = 1. Both for strong and 
weak drivings the system shows a normal Ohmic behavior, as 
it is depicted by the straight dashed line, which corresponds 
to ~ 5.2/ TV. The magnetization has an approximately linear 
profile, as shown in the inset, with a slight even-odd oscil- 
lation that is due to a staggered magnetic field. Data are 
obtained from MPO. 



lated in terms of the spin current, is observed and stud- 
ied also in the anisotropic Heisenberg spin chain. Our 
numerical data show that in the Ising-like regime, corre- 
sponding to anisotropy parameter |A| > 1, the system is 
a normal conductor for small driving fields, while it dis- 
plays negative differential conductivity for stronger driv- 
ings. The phenomenon is stable against the breaking of 
integrability, as, for instance, it persists in the presence 
of a staggered magnetic field. On the other hand, in the 
XY-like regime, |A| < 1, ballistic conduction with cur- 
rent proportional to driving field without negative differ- 
ential conductivity is observed at all drivings. Negative 
differential conductivity is schematically explained using 
the spectrum of one-magnon excitations, which become 
localized for |A| > 1 at the chain boundaries. 

The observed negative differential conductivity arises 
as an outcome of a beautiful interplay between co- 
herent many-body quantum dynamics of the inter- 
acting electrons/spin chain and incoherent charge 
injection-extraction/spin pumping operated by macro- 
scopic leads/spin baths. While our results are suggestive 
of a metal-insulator phase transition when driving our 
model systems far from equilibrium, new analytical ap- 
proaches are required to ascertain whether NDC survives 
at the thermodynamic limit. At any rate, the discussed 
magnon-localization mechanism for NDC is of potential 
interest for nanoscale devices such as current/heat diodes 
and transistor. 

Note added. After the completion of this work we be- 
came aware of a related paper— where it is shown that 
negative differential conductivity can be also achieved in 



a continuum limit of a fermionic tight binding model, 
driven far from equilibrium by embedding it between two 
free fermionic reservoirs. 
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APPENDIX A: NUMERICAL METHODS 

1. Quantum trajectories approach 

As compared to full density matrix simulations, the 
advantage of a quantum trajectory approach is that, in- 
stead of storing and evolving a density matrix of size 
TV x TV, one works with a stochastically evolving state 
vector of size TV (TV being the Hilbert space dimension 
of the system). The first two terms in the r.h.s. of 
Eq. (J?]) can be regarded as the evolution performed by 
an effective non-Hermitian Hamiltonian Tl e f{ = Tls + iJC, 
with JC = -jEm^mi the l &st term is responsi- 
ble for the so-called quantum jumps, as explained be- 
low. If the initial density matrix describes a pure state 
p(to) = |^(^o)) {0(^o) | •> after a small amount of time dt it 
evolves into the statistical mixture 



p(t + dt) = (l-J2 d P 



m I 90} \9o\ 



(Al) 

where dp m — (4>{tQ)\L\ n L m \(t)(tQ))dt and the new states 
are defined by 



-iTL c ¥vdt 



ffd * Wo)) 



V 1 - Em d Vr, 



i | fm 



I 077 



\L m \4>{t ))\\ 



(A2) 



Therefore, with probability dp m a jump to the state \(j} m ) 
occurs, while with probability 1 — ^ m dp m there are no 
jumps and the system evolves according to H c s- 

In practice, one starts from a pure random state |0(^o)) 
and, at intervals dt much smaller than the relevant dy- 
namical time scales, a random number e £ [0, 1] is cho- 
sen. If e < ^2 m dp m , the state of the system jumps to 
one of the states \<j) m ) (to \(f>\) if < e < dpi, to |</> 2 ) if 
dpi < e < dpi + dp2, and so on). On the other hand, 
if e > Y] m dp m the evolution with the non-Hermitian 
Hamiltonian 7Y c ff takes place, thus ending up in the state 
\4>o). This process has to be repeated as many times 
as n s top S = T/dt, where T is the total evolution time. 
Assuming that there exists a single out-of-equilibrium 
steady state p s , the expectation values (A) = tr (Ap s ) 
of any observable A are obtained after averaging in time 
(4>(t)\A\(j)(t)) up to a long enough time T, skipping the 
initial convergence time T conv needed for a solution of the 
Lindblad equation to converge into the stationary one. 
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In the cases discussed in this paper, the non-Hermitian 
evolution of e~*" eff< " has been simulated with a second 
order time Trotter expansion; a time slicing dt ~ 2 x 10 -3 
is sufficient to obtain accurate results. The fermionic 
systems taken into account are mapped into lattice spin 
models; simulations are then directly performed on spin- 
1/2 Hamiltonians [see Eqs. (fT4")l and (|18p]. We were able 
to simulate systems with sizes up to J\f ~ 10 5 (that is, 
N = log 2 N w 16 spin-1/2 particles). We calculated the 
stationary spin current (j) on the basis on Eq. Ijlip . That 
is, we computed (j) by summing up all down-up flips mi- 
nus all up-down flips at the right end of the chain or all 
up-down flips minus all down-up flips at the left end, and 
then dividing by the simulation time. In the fermionic 
language, this corresponds to counting the number of 
electrons that enter the system minus the number of elec- 
trons that leave the system per unit time through the 
right electrode (emitter) or the number of electrons that 
leave the system minus the number of electrons that enter 
the system per unit time through the left electrode (col- 
lector). We also checked that the current obtained in this 
way is, up to statistical fluctuations due to finite integra- 
tion times, equal to the one computed through Eqs. (|15[) . 
(fl9| . A good convergence for (j) is already reached at 
Tconv ~ 10 4 , while integration times T of one order of 
magnitude longer are required in order to determine the 
stationary magnetization profiles (e.g., T ~ 3 x 10 5 for 
N = 12 in Fig. El). 

Data presented for the Hubbard model and the spinless 
fermion model with N < 16 have been obtained with this 
approach. For further details of the implementation of 
the quantum trajectories approach see Ref. l47l 

2. Matrix Product Operator formalism 

We have used the MPO ansatz for spin-1/2 particles. 
An arbitrary density matrix p of a chain of N spins 1/2 
can be expanded over all possible products of Pauli oper- 
ators forming a basis of a 4 N dimensional Hilbert space 
of operators: 

|p>=5>-|a*>> (A3) 

£ 

where we used the compact notation a- = er^ 1 • • • a^ 1 , 
s = si ■ ■ ■ sat, and € {0, 1, 2, 3}, with <r° = 1, a 1 = 
a x ,(T 2 = <T y ,<T 3 — a z , and where lower indices in the 
Pauli operators denote the site number of the spin on 
which it operates. The use of ket notation in \p) outlines 
the fact that the density matrix p can be seen as a vec- 
tor in the operator Hilbert space spanned by the basis 
vectors \<J-). In the MPO ansatz, the expansion coeffi- 
cients Os are expressed as traces of products of N D x D 
dimensional matrices i = 1, . . . , N, as 

a.= tr(Af ■•■Aj*). (A4) 

Thus, a (density) operator (|A3[) is completely specified 
and thus parametrized in terms of a set of 4N matrices 



A*% four for each site i, Si = 0,1,2,3. We note that 
these matrices are not related to physically observable 
quantities. The propagator corresponding to the master 
equation ^ is written as a product of propagators for 
small steps of length dt, typically dt = 1CP 1 . Each small 
time-step propagator is then split using a third order 
Trotter expansion into parts composed of mutually com- 
muting two-spin terms. These nearest neighbor two-spin 
terms are then basic transformations performed within 
the MPO ansatz, namely after each such two-spin trans- 
formation, a singular value decomposition is performed 
in order to restore the shape of the ansatz (|A4[) . How- 
ever, this step has to be combined with a truncation of 
the resulting matrices to a smaller, fixed dimension D. 
Dimension D is then chosen as a parameter by which we 
control the accuracy of the method. Note that the mini- 
mal necessary dimension D is related to the bipartite en- 
tanglement of \p) in the Hilbert space of operators. This 
implies that MPO method will fail if the state p builds 
up strong quantum correlations over large distances. For 
a review on MPO techniques see Ref. [6l|, while details 
of the implementation of the MPO method in quantum 
master equations can be found in Ref. l33l 

Starting from an arbitrary initial density matrix p(to) 
we are interested in the asymptotic nonequilibrium 
steady state p s reached after a sufficiently long time 
of simulation, i.e., of relaxation. Once p s is obtained, 
various expectation values can be evaluated. In par- 
ticular, the spin current is obtained from Eq. Q19[) as 

(j) = Jx^[ps(^k a k+i ~ a k a k+i)]- The simulation time 
that is required to reach a stationary p s in the regime 
of NDC may greatly depend on the imposed gradient 
(driving field). In the most difficult situations that we 
studied here, the relaxation time after which the local 
current inhomogeneities decreased to a few percent was 
T ~ 4000. Note that this relaxation time is not directly 
comparable to the time needed in the quantum trajecto- 
ries approach, as there large T is needed also to perform 
statistical averaging. In the MPO approach averaging 
over Hilbert space is exact, up to the truncation imposed 
by finite matrix dimension D. The main advantage of 
MPO as compared to other approaches is that it gen- 
erally enables simulation of larger systems. For small 
drivings one can go up to N ~ 100, while in the present 
work, where systems are driven far from equilibrium, we 
could reach sizes N ~ 40. 

Data shown for the spinless fermion model with N > 
16 have been obtained with MPO technique. 



APPENDIX B: MAPPING FERMIONIC 
SYSTEMS INTO SPIN CHAIN MODELS 

In this appendix we sketch the main steps of the 
Jordan- Wigner Transformation (JWT)£2 mapping the 
Hubbard/i— V fermionic Hamiltonians with open bound- 
ary conditions into spin-1/2 ladder/chain models. In par- 
ticular, we show that the bath operators injecting or ex- 
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tracting fermions at the two borders of the system corre- 
spond, in the spin-1/2 picture, to operators nipping the 
edge spins. 

We start from the simpler case of the t — V model, 
described by Hamiltonian (fH))) . We first perform a JWT 
of fermions into hard-core bosons (aj, aj), defined by 



k=l 



(Bl) 



where rij — c^Cj is the fermion number operator. The 
operators aj for different sites commute, but they are 
not ordinary bosonic operators, since at most one bo- 
son is allowed on each site. One can indeed show that 



{a}) |0) = and, on the same site, {cij,aj} = 1. More- 
over, by using 



U(l-2n k ) JJ(l-2nfc/) = 1 



2rii 



(B2) 



k=l 



fc'=l 



[it follows from (1 — 2nk)(l — 2rik) — 1], and the fact 
that terms with different site index commute, we find 
rij = c-Cj = (Xjdj and c]c 3+ i = c](l — 2n,j)cj + i = oj-Oj+i. 

Then we have to transform the hard-core bosons into 
spin-1/2 particles, in a representation that identifies, at 
each site, the state |0) = a\l) with ||), and |1) = |0) 
with |t). If a" denote the corresponding spin-1/2 parti- 
cles (which of course obey the standard anticommutation 
rules {cr~ , at} — 6j } j>), the explicit mapping is given by: 



(B3) 



1. 



This leads quite straightforwardly to the following ex- 
pressions: 

t t t t 

c'jCj+i + c) +lCj = ajoj+i + a' j+1 aj 

Substituting them in Eq. (|16[) with open boundary con- 
ditions we finally get 

N-l N-l 



r 
7 



( 7 1 Z +^-(7V-l)-2^ ( 7j 



(B5) 



that is the Heisenberg Hamiltonian with J x = —t/2 and 
J z = V/4, plus an on-site transverse uniform magnetic 
field of strength V/2 and local transverse fields at the 
edge spins of the chain. 



The Lindblad operators are mapped into spin opera- 
tors that flip the outer spins. Indeed we have c\ = 
and ci = err on the left side; 



-N 



N 



(B6) 



and cm = e 17rN F ' aZ[ on the right side (where is the 
number of fermions in the leftmost N — k sites of the 
chain) . Since after action of the operator cr^ in (|B6|) the 
right-most site is always occupied, we can just as well 
rewrite Eq. (|B6p as 



TrJV" 



JV' 



(B7) 



where J\fp^ is the total number of spins up (occupied 



fermion states) in the entire lattice. 

Now, the crucial observation is that e l " J "p is an oper- 
ator which commutes (anticommutes) with the fcrmionic 
algebra consisting even (odd) number of fermionic op- 
erators. Since all the terms in our Lindblad equations 
conserve the parity of density operators^, i.e. they map, 
in the Liouville space sense, the products of even/odd 
number of fermionic operators Cj , cj into products of 
even/odd number of such operators, the Lindblad equa- 
tion ([4]) can be restricted to even-parity density operator 

subspace only. Thus, the factor e l7Tj ^F ' cancels from the 
Lindblad equation, even though the number of particles 
(or magnetization, in the spin language) is not conserved 
in an open system. Odd-parity density operator sub- 
space, which would include terms like |odd) (even|, etc., 
which change the number of fermions by an even num- 
ber, could be studied with a similar (but not identical) 
Lindblad equation with an additional minus sign in all 
the Lindblad terms. However, since all the physical ob- 
servables in question, say A = jk, rik, ■ ■ ., are represented 
as products of even number of fermionic operators, only 
the even-parity component of the density operator af- 
fects the expectation values (A) = tr [pA] . Thus we can 

safely forget the phase operator e l7rA ^ ' , and map Lind- 
blad operators from fermionic language to spin language 
in the Lindblad master equation (fj| as: — > &f and 
Cj — > aj , keeping in mind the issue of operator-space- 
parity in case expectation values of odd-parity operators 
would be needed. 

The charge current becomes then a spin current, while 
■i(l + (crj)) is the charge density on site j. We numerically 
checked that the addition of two local transverse fields on 
the border sites and of a uniform magnetic field does not 
qualitatively affect the NDC effect, as shown in Fig. [T7] 

We now show that the Hubbard model of Eq. (TT2")) can 
be mapped into the Heisenberg spin ladder of Eq. (|14p . 
Following the same steps for the t— V model, one first per- 
forms a double JWT of spin-up and spin-down fermions 
into two species of hard core bosons (aj, aj) and (&j, bj). 
Such transformation is defined by 



C 3;i 



(B8) 
(B9) 
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FIG. 17: (Color online). Spin current as a function of the 
driving strength for the Hamiltonian in Eq. (|B5j with t — 2, 
V = 8 (that corresponds to J x = 1, J z = 2, so that A = 2). 
In contrast to the standard XXZ model, here we added a 
uniform transverse magnetic field of intensity V/2, plus two 
transverse fields of strength — V/4 at the border sites of the 
Heisenberg chain. The system-bath coupling is set equal to 
r = 4; the simulation time (QT approach) is T = 2.5 x 10 4 . 



Then the hard core bosons are transformed into two 
species of spin- 1/2 particles (<7j and r. 



jj, in a repre- 
sentation analogous to that of Eqs. (|B3|) . Proceeding 
along the same transformations as before, and using 

n j,1 n j,l — 5 (°j + 1) | ( T j + l)i we nna lly arrive at the 
spin ladder Hamiltonian (|14[) . In complete analogy with 
the t — V model, the Lindblad operators of Eqs. (O, ([6]) 
are mapped into spin operators flipping the outer spins 
of the ladder. 



APPENDIX C: SPIN-SPIN CORRELATIONS 

We provide here some numerical data on the behavior 
of the steady-state spin-spin correlation functions for the 
Heisenberg chain driven far from equilibrium, aimed at 
investigating the emergence, in the gapped regime and 
at strong driving strength /, of a long-range correlation 
order—. The steady-state spin-spin correlation function 
is defined as: 



(CI) 



where (ofcrj) = tr (erf 0jp s ) and (erf) = tr (crfp s ), while 
averages are taken on the steady state. Results are plot- 
ted in Fig. |TH] With increasing / and above /*, a dra- 
matic slowing down of the correlation decay is appar- 
ent. Unfortunately, a more quantitative understanding 
of a possible critical behavior in the large-/ regime is at 
present out of our capabilities, since it would require the 
analysis of much bigger systems, and this is numerically 
not accessible for the model under consideration. 
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FIG. 18: (Color online). Spin-spin correlation function C(i, j) 
in a chain with 40 spins, for / = 0.1 (top row) and / = 0.9 
(bottom row); here we set A = 2 and T = 4. The code in 
the left panels denotes \og 10 [C(i,j)]. The right plots display 
correlation function C(r = \i — j\) along the diagonal denoted 
by a dashed line in the left plots; the case with N — 20 spins 
is also shown. Data are obtained using the MPO ansatz. 



APPENDIX D: SOLUTION OF THE 
ONE-MAGNON MODEL 

The matrix associated with the Hamiltonian (|18p in 
the one-magnon basis {|1) , |2) , | N)} is tridiagonal and 
reads as follows: 



H 



r a' (3 

a (3 
(3 a (3 



(3 a (3 
p a' 



(Dl) 



where a' = (N-3)J Z , a = (N-5)J Z , and = 2J X . Note 
that the overall magnetization is conserved by the XXZ 
Hamiltonian, so the states corresponding to the other 
spin sectors are not coupled (by the autonomous evolu- 
tion) to the one-magnon sector. The eigenvalues of H 
are given by the roots of the characteristic polynomial 
D N {E) = det(El - H). 

First, it is convenient to solve the eigenvalue problem 



(0) 

N 



D 



N 



The following 



for a — a'. We define D 
recurrence relation holds 

Df =(E- a)Df\ - p 2 Df} 2 . (D2) 

The general solution to this difference equation is 

D<W = {3 m (Ae lm0 + Be-™*) , (D3) 

with 



cos ( 



E 



2j3 



(D4) 
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and the constants A and B determined from the condi- 
tions 



D 



We finally obtain 



(0) 



(o) 



(E-a)D^> -/3 2 . 



D 



N 



whose solutions are 



E$ = a + 2/3 cos 



(0) _ P N sin[(/V + l)0] 



(D5) 



(D6) 



JV + 1 



, m = l,...,N. (D7) 



These eigenvalues are located in the energy band a — 
2\(3\ < E < a + 2\P\. Therefore, <f> is always real and the 
eigenstates ( "Bloch orbitals" ) 



(0)\ 



nmk . 

sin | ^r-^ 1 , (D8) 



JV + 1 ^ 1 

k— 1 



corresponding to the eigenvalues E$ are delocalized 
along the spin chain. 

In the case a ^ a' we obtain 



D N = Df +2{a- ( J)Df_ l + {a-a') 2 Df_ 2 



N 



{sin[(iV + 1) 



(D9) 



2Asin(JV</>) + A 2 sin[(JV - !)</>}} , 



where we have used 



a' - a _ Jz_ _ . 
~ J, ~ ' 



(D10) 



The equation — can be analytically solved for large 
N. There always exist at least N — 2 delocalized solutions 
lying in the energy band between a — 2\/3\ and a + 2\f3\. 
The "molecular orbitals" \ip±) w -j^(\ip L ) ±\ip R )) appear 
when |A| > 1. 

If A > 1, we have </> = i\ (x G K), e x w A, 



1 - e- 2 x 
1 - e- 2Ar * 

1 - e- 2 ^ 
1 - e" 27V x 



N 

E« 

m=l 
N 



-(m-l)x 



-(JV-m)x 



to) . (Dll) 



The states I^L.fl) ar e centered at sites 1 and iV, respec- 
tively, and their localization length £ w l/\ = l/ln(A). 
The corresponding eigenvalues are given by 



Ei « E N w a + /3 A 



(D12) 



E-a 
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FIG. 19: (Color online). The one-magnon model for AT = 
16 spins. Left panel: energy spectrum as a function of A; 
the large- N analytic result (|D12|I is shown by thick dashed 
lines. Right panel: molecular states \ip±) = X^m=i( Cm )± l m )> 
at A = 2; the numerically computed coefficients (c m )+ and 
(cm)_ are shown as circles and triangles, the large- iV analytic 
results as full and dashed curves. 



If A < -1, <j> = ix + 7T (x G K), ex w -A, 



A' 



m— 1 

1 _ _2 X W 



1 - e- 2Ar x 



These states have localization length £ s» 1/x = 
l/ln(— A). The corresponding eigenvalues are again 
given by Eq. l(Dl2"]) . 

The gap 6e between the energy levels E\ and En 
shrinks exponentially with the system size: 



N 

Se oc \(tp L \ipR)\ ~ exp ( — - 



(D14) 



so that the coherent tunneling between sites 1 and JV re- 
quires a time scale which grows exponentially with JV. 
Therefore, for the purposes of our present investigation 
we can say that a spin-flip excitation created at one 
boundary of a ferromagnetic domain remains in practice 
exponentially localized over a length i = 1/ In |A|. 

Numerical illustrations, for a chain of JV = 16 spins, 
of the appearance, for A| > 1, of two states outside the 
energy band (a — 2\(3\,a + 2\f3\) and of the molecular 
orbitals \ip±), are provided in Fig. fT9l 

Finally, we note that the one-magnon localization 
model discussed in this section has deep similarities with 
a tight-binding model discussed in the context of surface 
physics^. 
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